home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / recurse.h < prev    next >
Encoding:
C/C++ Source or Header  |  2001-11-01  |  7.7 KB  |  126 lines
  1. /* specfunc/recurse.h
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author: G. Jungman */
  21.  
  22. #ifndef _RECURSE_H_
  23. #define _RECURSE_H_
  24.  
  25. #define CONCAT(a,b) a ## _ ## b
  26.  
  27.  
  28. /* n_max >= n_min + 2
  29.  * f[n+1] + a[n] f[n] + b[n] f[n-1] = 0
  30.  *
  31.  * Trivial forward recurrence.
  32.  */
  33. #define GEN_RECURSE_FORWARD_SIMPLE(func)                                      \
  34. int CONCAT(recurse_forward_simple, func) (                                    \
  35.                                const int n_max, const int n_min,              \
  36.                                const double parameters[],                     \
  37.                                const double f_n_min,                          \
  38.                    const double f_n_min_p1,                       \
  39.                                double * f,                                    \
  40.                    double * f_n_max                               \
  41.                                )                                              \
  42. {                                                                             \
  43.   int n;                                                                      \
  44.                                                                               \
  45.   if(f == 0) {                                                                \
  46.     double f2 = f_n_min;                                                      \
  47.     double f1 = f_n_min_p1;                                                   \
  48.     double f0;                                                                \
  49.     for(n=n_min+2; n<=n_max; n++) {                                           \
  50.       f0 = -REC_COEFF_A(n-1,parameters) * f1 - REC_COEFF_B(n-1, parameters) * f2; \
  51.       f2 = f1;                                                                \
  52.       f1 = f0;                                                                \
  53.     }                                                                         \
  54.     *f_n_max = f0;                                                            \
  55.   }                                                                           \
  56.   else {                                                                      \
  57.     f[n_min]     = f_n_min;                                                   \
  58.     f[n_min + 1] = f_n_min_p1;                                                \
  59.     for(n=n_min+2; n<=n_max; n++) {                                           \
  60.       f[n] = -REC_COEFF_A(n-1,parameters) * f[n-1] - REC_COEFF_B(n-1, parameters) * f[n-2]; \
  61.     }                                                                         \
  62.     *f_n_max = f[n_max];                                                      \
  63.   }                                                                           \
  64.                                                                               \
  65.   return GSL_SUCCESS;                                                         \
  66. }                                                                             \
  67.  
  68.  
  69. /* n_start >= n_max >= n_min 
  70.  * f[n+1] + a[n] f[n] + b[n] f[n-1] = 0
  71.  *
  72.  * Generate the minimal solution of the above recursion relation,
  73.  * with the simplest form of the normalization condition, f[n_min] given.
  74.  * [Gautschi, SIAM Rev. 9, 24 (1967); (3.9) with s[n]=0]
  75.  */
  76. #define GEN_RECURSE_BACKWARD_MINIMAL_SIMPLE(func)                             \
  77. int CONCAT(recurse_backward_minimal_simple, func) (                           \
  78.                            const int n_start,                             \
  79.                                const int n_max, const int n_min,              \
  80.                                const double parameters[],                     \
  81.                                const double f_n_min,                          \
  82.                                double * f,                                    \
  83.                    double * f_n_max                               \
  84.                                )                                              \
  85. {                                                                             \
  86.   int n;                                                                      \
  87.   double r_n = 0.;                                                            \
  88.   double r_nm1;                                                               \
  89.   double ratio;                                                               \
  90.                                                                               \
  91.   for(n=n_start; n > n_max; n--) {                                            \
  92.     r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
  93.     r_n = r_nm1;                                                              \
  94.   }                                                                           \
  95.                                                                               \
  96.   if(f != 0) {                                                                \
  97.     f[n_max] = 10.*DBL_MIN;                                                      \
  98.     for(n=n_max; n > n_min; n--) {                                               \
  99.       r_nm1  = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
  100.       f[n-1] = f[n] / r_nm1;                                                     \
  101.       r_n = r_nm1;                                                               \
  102.     }                                                                         \
  103.     ratio = f_n_min / f[n_min];                                               \
  104.     for(n=n_min; n<=n_max; n++) {                                             \
  105.       f[n] *= ratio;                                                          \
  106.     }                                                                         \
  107.   }                                                                           \
  108.   else {                                                                      \
  109.     double f_nm1;                                                             \
  110.     double f_n = 10.*DBL_MIN;                                                 \
  111.     *f_n_max = f_n;                                                           \
  112.     for(n=n_max; n > n_min; n--) {                                               \
  113.       r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n);  \
  114.       f_nm1 = f_n / r_nm1;                                                       \
  115.       r_n = r_nm1;                                                               \
  116.     }                                                                         \
  117.     ratio = f_n_min / f_nm1;                                                  \
  118.     *f_n_max *= ratio;                                                        \
  119.   }                                                                           \
  120.                                                                               \
  121.   return GSL_SUCCESS;                                                         \
  122. }                                                                             \
  123.  
  124.  
  125. #endif /* !_RECURSE_H_ */
  126.